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Abstract 


This paper presents a reliability- and robustness-based formulation for robust control 
synthesis for systems with probabilistic uncertainty. In a reliability-based formula- 
tion, the probability of violating design requirements prescribed by inequality con- 
straints is minimized. In a robustness-based formulation, a metric which measures 
the tendency of a random variable/process to cluster close to a target scalar /function 
is minimized. A multi-objective optimization procedure, which combines stability 
and performance requirements in time and frequency domains, is used to search 
for robustly optimal compensators. Some of the fundamental differences between 
the proposed strategy and conventional robust control methods are: (i) unneces- 
sary conservatism is eliminated since there is not need for convex supports, (ii) the 
most likely plants are favored during synthesis allowing for probabilistic robust op- 
timality, (iii) the tradeoff between robust stability and robust performance can be 
explored numerically, (iv) the uncertainty set is closely related to parameters with 
clear physical meaning, and (v) compensators with improved robust characteristics 
for a given control structure can be synthesized. 

Several numerical methods for estimation, including the Hammersley sequence sam- 
pling method, the First Order Reliability method, and the First- and Second- 
Moment-Second-Order-Methods, are compared. Examples using output-feedback 
and full-state feedback with state estimation are used to demonstrate and validate 
the methodology. 
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1 Introduction 


Achieving balance between stability and performance in the presence of uncertainties 
is one of the fundamental challenges faced by control engineers. Tradeoffs must be 
made to reach acceptable levels of stability and performance with adequate robust- 
ness to parameter uncertainty. These tradeoffs are explicitly linked to the control 
engineer’s choice of uncertainty model as well as how that model is exploited in the 
synthesis process. Usually, the assumed uncertainty model has a profound impact 
on the performance robustness of the closed-loop system. 

Several uncertainty models, such as norm-bounded perturbations, interval anal- 
ysis, fuzzy sets and probabilistic methods [1-3] are typically used. The most com- 
monly used robust control methods [4] are /i-synthesis and iUinfinity. In these 
methods, uncertainty is modeled with norm-bounded complex perturbations of ar- 
bitrary structure about a nominal plant. This treatment is used primarily because 
it leads to a tractable set of sufficient conditions for robust stability, making the ap- 
proach computationally efficient. These methods are based on the most pessimistic 
value of performance among the possible ones, usually referred to as “worst-case”. 
This worst-case performance is usually realized only by a single member of the 
uncertain model set and by a particular input signal. No information is provided 
regarding the likelihood that this worst-case will ever occur in practice. In addition, 
the intrinsic mathematical requirements of the approach usually lead to conservative 
models of uncertainty, over-conservative designs and complicated compensators. 

Probabilistic uncertainty not only defines a set of plants where the actual dy- 
namic system is assumed to reside but also associates a weight, the value of the 
probability density function, to each member of the set. In contrast to conventional 
robust control methods, this “additional dimension” allows the pursuit of robustly 
optimal solutions in the probabilistic sense. For instance, reliability-based design 
searches for solutions that minimize the probability of violating design requirements 
prescribed in terms of inequality constraints. Hence, reliability-based control design 
searches for the compensator that places as much probability as possible within 
the region where the design requirements are satisfied. Notice that this allows the 
search for the compensator with the best robustness for a given control structure, 
e.g., the most robust PID controller, even though the violation of some the design 
requirements for some of the plants in the uncertainty set is possibly unavoidable. 

Synthesis approaches based on random searches [5-7] and stochastic gradient al- 
gorithms [8-10] have been applied to probabilistic robust control. In these studies, 
random sampling is the primary tool for assessing and pursuing acceptable levels 
of robustness in the control solution. On the other hand, asymptotic approxima- 
tions [11, 12] for the estimation of failure probabilities have only been used as a 
control analysis tool. The works [5-7, 13, 14] and references cited therein are espe- 
cially relevant to this paper. Even though they lay down the basic framework for 
the reliability control synthesis of engineering problems, important aspects of the 
formulation and of the solution method remain to be explored and refined. This 
article addresses and extends some of those aspects. 
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The main contributions of this study to the state of the art in the subject are 
as follows: 

1. The use of robustness-based metrics for minimizing the performance degrada- 
tion caused by uncertainty. 

2. The use of bounds on the reliability metrics based on the first two order 
moments. This practice reduces considerably the computational demands of 
the synthesis algorithms based on the estimation of failure probabilities. 

3. The use of shapable failure domains within the reliability formulation. This 
allows the integration of robustness considerations into the conventional reli- 
ability approach. 

4. The integrated use of deterministic sampling and asymptotic approximations 
for the estimation of reliability metrics. This hybrid approach (i) reduces the 
computational complexity of the synthesis algorithm without compromising 
the accuracy of the results, (ii) eliminates the random character of the esti- 
mation, and (iii) eliminates the high computational demands associated with 
the estimation of small failure probabilities via sampling. 

This paper is organized as follows. Section 2 presents basic concepts related to 
control and probabilistic uncertainty. Section 3 introduces reliability metrics for 
random variables and processes and presents realizations to stability and time- and 
frequency-dependent performance metrics. Mean and variance based bounds to the 
reliability metrics are also derived therein. Robustness-based metrics for random 
variables and processes are introduced in Section 4. Section 5 presents the numer- 
ical methods used to estimate the above mentioned metrics. The control synthesis 
procedure is presented in Section 6, where specifics of both the reliability and the 
robustness-based formulations are examined. Two examples are presented in Sec- 
tion 7, where a satellite’s attitude control problem and the disturbance rejection in 
a flexible beam are used to demonstrate the method. Finally, some conclusions are 
stated in Section 8. 

2 System Dynamics 

Let p be a vector of random variables used to model the uncertain parameters of 
the system. In this study, p is prescribed a priori by the joint probability density 
function (PDF) / p (p) or equivalently by the cumulative distribution function (CDF) 
Fp(p) 1 . The set of values that p could take, called the support of p, will be denoted 
as A p . 

Consider the probabilistic model A4(p) of a Linear Time Invariant (LTI) system, 
where the dependence of the model on the uncertain parameters could be non-linear. 
The reader must notice however, that the developments presented herein do not 
require the system to be LTI. The propagation of A p through A4 leads to a set 

1 In these expressions, the subscript refers to the symbol used for the random variable while the 
value in parenthesis refers to a particular realization of it. 
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of uncertain plant models in which the physical system is assumed to reside. The 
probability distribution of a plant within this set is fully determined by Al(p) and 
/p(p). In a transfer function representation, we will refer to the uncertain plant as 
G'(p) and to the compensator as K( k), where k is the vector of design parameters 
to be determined. Alternatively, a state space realization of A! (p) leads to 

x = A(p)x + B(p)u + F(p)z (1) 

y = C(p)x + D(p)u + E(p)v (2) 

where x is the state, u is the control, z is process noise, y is the system output and v 
is sensor noise. The noise signals are commonly modeled as delta correlated Gaussian 
white noises satisfying E[z] = 0 and E[z(t)z T (t + k)] = S<5(ac) , where z = [z T , v 7 ] T , 
S is the spectral density matrix and E [•] is the expected value operator. In what 
follows, the explicit dependence on p is omitted while D is assumed to be zero. 

Important properties used in control design, such as pole placement and the 
Separation Principle, do not hold due to the offset between the deterministic math- 
ematical model and the actual dynamic system. The effects of parametric uncer- 
tainty on the Separation Principle are considered next. For the full-state feedback 
law u = — Gx and a full-order observer with gain L based on the expected plant 
E[Af(p)] (any other deterministic plant such as ,A/f(E[p]) could be used instead), 
the observer equation and the closed-loop dynamics for velocity feedback are given 

by 


X = (E[A] - E[B]G + L(C - E[C])) x 

(3) 

x = Ax + Bz 

(4) 

y = Cx + Ez 

(5) 


A - BG 

BG 

A — E[A] + (E[B] - B)G 
+L(E[C] - C) 

(B — E[B])G + E[A] — LE[C] 


F 

0 

0 

-LE 


where x is the estimation of x, x = [x r , g 7 ] T is the augmented state vector, g = x— x 
is the estimation error, C = [C 7 |0 7 ] T and E = [0 r |E 7 ] 7 . The vector k is formed by 
the feedback gain G and the observer gain L. Notice that the separation principle 
holds, i.e., A is upper triangular, if the deterministic plant used to generate the 
observer matches the dynamic system exactly. However, in general, the Separation 
Principle does not hold due to uncertainty in the plant model. In addition, the 
random closed-loop poles do not occur at the locations selected for the full-state 
feedback, i.e. poles of the Api subsystem, nor at the locations for the full-order 
observer, i.e. poles of the A2,2 subsystem. 
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In this paper we propose a methodology for robust control design for systems 
with probabilistic uncertainty. Some of the major differences between the proposed 
strategy and conventional robust control methods are: (i) unnecessary conservatism 
is eliminated since there is not need for convex or bounded supports, (ii) the most 
likely plants are favored during synthesis allowing for probabilistic robust optimality, 

(iii) the tradeoff between robust stability and robust performance can be explored, 

(iv) the uncertainty set is closely related to parameters with clear physical mean- 
ing, and (v) compensators with improved robust characteristics for a given control 
structure, e.g., PID, can be designed. 


3 Reliability-Based Metrics 

The propagation of a fixed set of parameters of the plant through conventional con- 
trol analysis tools leads to set of scalar quantities, e.g., closed loop poles, and a 
set of functions, e.g., step responses and Bode plots. The propagation of proba- 
bilistic uncertainty through the same tools leads to random variables, e.g., random 
closed- loop poles, and random processes, e.g., the step responses become random 
processes parameterized by time and the Bode plots become random processes pa- 
rameterized by frequency. In this section we first introduce reliability metrics for 
random variables and processes. Such metrics will be used to quantify the violation 
of the design requirements. Specific realizations corresponding to stability, time 
and frequency requirements are then provided. In general, we will use x and x(h) to 
denote a random variable and a random process dependent on p through the plant 
model. For the random process x(h), h refers to an arbitrary variable such as time 
or frequency. 

3.1 Random Variables 

We first introduce the concept of probability of failure. Let x(p) be the random 
variable of interest. Let x > x be a design requirement. The event x < x will 
be referred to as failure. The corresponding failure set is given by T = {x | x £ 
(— oo,x]}, where the failure boundary x is a deterministic quantity prescribed in 
advance. The admissible domain, namely A = {x \ x G (x, oo)}, is the complement 
of the failure domain. The same type of discrimination can be done in the parameter 
space p by using x(p). The function g(p,x) = x(p) — x, called the limit state 
function, divides the parameter space in two parts, the domain leading to A, i.e., 
g{ p,x) > 0, and the domain leading to J~, i.e., g{ p,x) < 0. Hence, T results from 
mapping the set {p € A p | g(p,x) < 0} through x(p). In this case, the probability 
of failure Pf is given by 


p f = p[ x < x] = f f x (0d£ = f f P (p)dp (6) 

J£<x Jg < 0 

Similar expressions can be derived if the design requirement is x < x. These expres- 
sions describe reliability metrics for the random variable x when a single constraint 
is present, i.e., x > x or x < x. A reliability metric for the random variable x having 
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both constraints can be easily formed 


r x {x,x) = r x (x) +f x (x) 

(7) 

where 


r x (x) = P[x < x\ = F x (x) 

(8) 

f x (x) = P[x > x\ = 1 - F x (x) 

(9) 


Notice that r x (x) is equivalent to Pf in Equation (6). We will refer to x and x 
as the boundaries of the failure domain F = {x \ x E (— oo,x] U (x, oo)}. Notice 
that the under-bar and the over-bar refer to the bound from below and the bound 
from above of the admissible domain A = {x \ x E \x,x)}. This convention will be 
used for the remainder of the paper. Notice that the mapping of the corresponding 
limit state function through x(p) leads to the failure boundary(s). Hence, there is 
a direct correspondence between T and g. A sketch with relevant information is 
provided in Figure 1. 



Figure 1. Sketch of the reliability metric for x. 


3.2 Random Processes 

The random process x(h) can be considered as the parameterization of the random 
variable x by the deterministic variable h. In this paper h E [0, oo] is assumed. The 
random process x(h) is specified by the set of CDFs [15] F x ^(x,h). For instance, 
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the system output y{t) is prescribed by F y ^ t \(y,t). The evaluation of the process at 
a particular h value, say hi , leads to the random variable x{h{) whose CDF is given 
by F x (x) = F x (h)(x, hi). In general, the support and the percentiles of x(h) vary 
with h. 

Let x(h) > x(h) for h € and x(h) < x(h) for h € [^ 3 ,^ 4 ] be design 

requirements for the random process x(h). In this paper, reliability metrics for 
processes are formulated by extending the ideas presented above. This is attained 
by integrating the reliability metric in Equation (7) for the random variable x(hi) 
in the h— interval of interest. In this context, a reliability metric for x(h) is cast as 

r x (h) (*(■)>*(■)) = !*(.)(*(■)) J rr x{h) (x(h)) ( 10 ) 

where 

£*(/») (*(■)) = “ J h p i x ( h ) < x(h)]dh = — l _ — F x{h) (x(h), h)dh (11) 

r x (h)( x ( 0 ) = h ~ ~ J P [x{h) > x(h)]dh = - - — — J 1 - F x ^(x(h),h)dh 

(12) 

are the costs of violating the lower and upper constraints respectively. These con- 
straints, namely x(h) and x(h), will also be referred to as failure boundary functions. 
Notice that the failure domain 


F = l (J {(x,h)\x < x(h)} > U < |J {(x,h)\x>x(h)} 

J ^ /ie [ft .3 ,/i4] 

is delimited by the failure boundaries. Observe that Equation (10) is a natural 
extension of Equation (7). A sketch with some of the pertinent metrics is provided 
in Figure 2. On the top plot, the 1, 25, 75 and 99 percentiles 2 are shown along with 
the failure boundaries x(h) and x(h) which, for this example, are linear. On the 
bottom plot, the integrands of Equations (11-12) corresponding to the configuration 
in the top plot are shown. Notice that if the process is contained within the set A, the 
reliability metric r x n^ is zero, meaning that the inequality constraints are satisfied 
for all parameter values in A p . 


3.3 Realizations 


3.3.1 Robust Stability 


A LTI system is robustly stable if all its poles are in the open left half of the complex 
plane for all possible values of the random parameters. A reliability assessment of 
stability is given by 


P 


IM*]>°) 

,i=i 


e 


^Recall that the m percentile, given by the x values satisfying F x ^{x,h) = m/100, defines a 
line under which m% of the probability lies. These lines allow us the visualize the h dependence of 
the PDF. 
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X 



99 percentile 
J- 75 percentile 
... 25 percentile 
1 percentile 


h 3 h/ 



Figure 2. Sketch of the reliability metric for x(h). 


where Si with i = 1, 2, . . . v is a random pole, 3i[-] is the real part operator and e is 
the resulting probability of instability. Robust stability is attained if e = 0. Stability 
can also be cast via 

A = max{5i[si] , 3? [ 52 ] , • • ■ , 3?[s„] } (13) 

In terms of A, the probability of instability is given by r\( 0). Robust stability 
is attained if r,\(0) = 0. Several comments are now pertinent. Reaching robust 
stability may not be feasible for the given support A p (even though it is bounded) 
and the assumed control structure R(k). Notice also that the acceptance of a 
small non-zero probability of instability could be desirable from the performance 
point of view. For instance, by allowing the right low-probability tail of f\ (A) to 
lie on the closed right half of the complex plane, a significant enhancement in the 
performance of the plants associated with the high probability portions of the PDF 
can be attained. Rather than advocating for the acceptance of the risk that this 
practice implies, we would like to highlight that the tradeoff between robustness and 
performance can be studied by allowing small positive values of e. 

3.3.2 Time-Domain 

Quite frequently, performance requirements are prescribed in terms of time-domain 
specifications. The propagation of / P (p) through the system dynamics leads to 
random processes for the time responses. Denote by x[t) an arbitrary random 
process with CDF F x u\(x,t). Such process is parameterized by p, time t and the 
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compensator design variable k. The dependence of x(t) on k has been omitted 
for the sake of simplifying the notation. Reliability metrics for relevant processes 
can be cast using Equation (10). For instance, while settling time and overshoot 
requirements for a Single Input Single Output (SISO) system, i.e., u e M and y G M, 
are integrated using r y ^(y(t),y(t)), the control saturation requirement |u| < u max 
leads to Umax , U rn ax ) • 

A reliability metric for assessing the effects of noise on the uncertain plant is 
formulated next. The state covariance matrix, defined as Q(t) = E[x(i)x T (i)], is 
given by the solution to the Lyapunov equation 

Q = AQ + QA r + BSB T (14) 

subject to Q(0) = Qo- The output covariance, defined as E[y(t)y i (f)], reaches the 
steady-state Root Mean Square (RMS) value 

r~ ~ 1 1/2 

y rms = lim diag CQ(f)C r (15) 

t — >oo L 

Notice that uncertainty in p makes y rms a random vector. For instance, a relia- 
bility metric that penalizes the violation y rm s > y rms where y rm s £ M is given by 

rms (y rms)' 

3.3.3 Frequency-Domain 

The propagation of / P (p) through the system dynamics onto the frequency domain 
leads to random processes of the form x(uj), whose probabilistic behavior at each 
lo is specified by F x r^{x,uj). Here, x{uj) is any real frequency dependent metric 
of the feedback loop, e.g. Bode magnitude. This random process is parameterized 
by p, frequency to and the compensator design variable k. A reliability metric for 
x(u)) is r x r U) )(x( to),x(u)). For instance, conventional control requirements [16] for 
disturbance rejection, noise attenuation and reference tracking can be cast in terms 
of the loop transfer function q( to) = |GAT|. Low frequency requirements can be cast 
using r q ^(q(to)) with q(uj) = 1 and high frequency requirements with r q ^(q(io)) 
for which q(uj) has a proper roll off, each over a suitable range of frequencies. 

3.4 Reliability Bounds 

The following lemma, based on the Tchebyshev inequality, sets bounds for the reli- 
ability metrics introduced above in terms of the first and second order moments. 

Lemma 1. Let E[-] and V[-] denote the expected value and variance operators. 

-*(-) - (E[^ -a) 2 if ~ < (16) 

V[x\ 

r x (x) < -j— ymw */* > E N 

(a — E[aJ) z 

8 


(17) 



Proof. If e > 0 

V[s] > [ (£ - E[x]) 2 f x (€)d£ > e 2 P [\x - E[s]| > e] > e 2 P [x < E[x] - e] 

^|C-E[*]|> e 

We obtain Equation (16) using x = E[,x] — e. Equation (17) is derived following the 
same lines. □ 


Notice that these bounds apply to any PDF of x. Since the bounds are solely 
dependent on the first two order moments they can be estimated more accurately 
than exact failure probabilities. Notice however that the use of bounds instead of 
reliability metrics introduces conservatism into the solution. A metric, based on 
these bounds, is defined as 


b x (x 

,x) =b x {x 

) + b x (x) 

(18) 

where 

A 1 

r vw 

if x < E[x] 
otherwise 


h(x) = | 

1 (E[x]-x ) 2 

1 OO 

(19) 

b x (x) = | 

r vm 

if x > E[x] 
otherwise 


1 (x-E[x ]) 2 

1 OO 

(20) 

As before, the under-bar and over-bar on b refer to the way 
event is defined. Notice that r x < b x , r x (x) < b x (x) and f x (x ] 

in which the failure 
) < b x (x). The same 


idea, extended to random processes, leads to 


K(h) (*(■)>*(■)) = K(h ) (*(•)) + b x (h)(x{-)) 


(21) 


where 


hx 4 J h^Ki In l K(x{h))dh if x(h) < E[s(/i)] V/i G [h u h 2 \ ^ 

~ x( ' ~ I oo otherwise 


b x(h] (x{-)) = J SWSs fha b x(x{h))dh if x{h) > E[x{h)} Mh € [h 3 , h 4 ] ^ 

x 1 ^ otherwise 


oo 


As before, r x{h) < b x{h) , < b x(h) (x{-)) and r x{h) (x(-)) < b x(h )(x(-)). 


4 Robustness-Based Metrics 

Some performance requirements might not be properly captured in a conventional 
reliability formulation since they focus on the bulk portion of the PDF. Throughout 
this paper the term robustness is used to describe the design characteristic that mea- 
sures the performance degradation from an ideal deterministic behavior caused by 
uncertainty. Robustness metrics, that quantify such a characteristic, are presented 
next. For random variables, the index 

r x (x) = [ (£ - x) 2 f x (t)d£, = V[x] + E[x] [E[x] - 2x) + x 2 (24) 

J A x 
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is a measure of the concentration of f x {x) about the deterministic target value x. In 
the ideal case when f x (x) = 5(x — x), we obtain r x (x) = 0. For random processes, 
the index 

T x(h) ($(■)) = /?6 ]_ J h T x(x{h))dh (25) 

is a measure of the concentration of the process x(h) about the deterministic target 
function x(h) for h £ [h§,ho[. Notice that the evaluation of the above expressions 
only requires of the first two order moments of the process. 

Robustness-based metrics parallel to the ones provided in Section 3.3 can easily 
be posed. For instance, if y rm s is the RMS steady state value of an error signal, the 
index T yrms (0) quantifies the offset between the target behavior y rms = 0 and the 
random variable y rms . Likewise, the metric (0) quantifies the offset between the 
random process u (t) £ M and the target u = 0, for which no actuation is required. 

5 Numerical Estimation 

Methods used for the estimation of the reliability- and robustness-based metrics in- 
troduced above are presented herein. Historically, mathematical studies of reliability 
engineering systems have focused on approximation of “probability of failure.” Some 
of the techniques historically used for this approximation are particularly suited for 
estimating probabilities near zero, which is the range in which one hopes a prob- 
ability of failure will lie. In this study, we will be concerned with estimating the 
probabilities of various random events. In this paper, events in a reliability frame- 
work might be referred to as “failure events” to be consistent with previous usage 
whether they actually represent a failure of some sort or not. The reader should also 
be aware that any estimation technique that works well for failure probabilities near 
zero also work well for events of probability near one after using the complementary 
event instead. 

Only the estimation of statistics for x will be addressed since the same tools can 
be extended to processes. Such an extension is as follows. For the random process 
x(h), create a uniform sample of e points hi, / 12 , • • ., h e in the h-domain. We will 
refer to these samples as the e h-samples. Statistics for the resulting e random 
variables Xi = x(hi), i = 1,2 ... e are then used to estimate the pertinent integrals, 
i.e. , Equations (11), (12), (22), (23), and (25). 

5.1 Hammersley-Sequence-Sampling (HSS) 

HSS generates representative deterministic samples of / p (p). The error of approxi- 
mating an integral by a finite number of samples of the integrand, depends on the 
uniformity of the points used to generate the samples rather than on their random- 
ness. This fact has motivated the development of deterministic sampling techniques 
such as HSS. While conventional Monte Carlo Sampling (MCS) is based on the gen- 
eration of random points on the unit hypercube, HSS is based on the generation of 
an evenly distributed set of points. The n Hammersley samples are generated by 
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transforming the n Hammersley points m* with i = 1, 2, ... n through the inverse 
CDF of the uncertain parameter 

Pi = fp'K) (26) 

The Hammersley points [17] are generated as follows. Let R > 1 be an integer. The 
integer i expressed in radix- R notation is given by 

V 

* = E i i RJ 

3=0 

where v = ffoorjlog^i} and ij is an integer between 0 and R — 1 for each j = 
0, 1, . . . v. The inverse radix number of i, namely <f>R(i), is given by 

M*) = E ih 

3=0 

The Hammersley points in a k-dimensional unit hypercube are given by the sequence 

m, = 1 - [*/ n, (fiRi (*) > 4>R2 (*)>••• ^Rk-i (i)\ T ( 27 ) 

where i = 1, 2 . . . n and Rj is the jth prime number. 

HSS requires far fewer samples [18] than MCS [5-7, 14] for a given confidence 
level. Improvements by a factor of three to one hundred in the convergence rate 
of the estimated first two order moments have been reported [19]. Another way to 
observe the advantages of HSS over MCS is to note that, for the same sample size, 
the HSS sample is more representative of the random distribution than the MCS 
sample. Figure 3 shows a comparison between HSS and MCS. At the top, n = 200 
points on the unit hypercube are shown. At the bottom, the corresponding samples 
for / P (p) = f a (a)fb(b), where f a (a) = N(0,1 ) and f b (b) = 5(3,2) with A b = [0,1] 
are displayed. Here, N and B denote a Normal and a Beta distribution. Substantial 
differences in the uniformity of the points and in the clustering of the samples are 
observed. In addition, the results of calculations based on HSS samples are not 
random. On the other hand, numerical experiments performed using MCS show 
variability due to such factors as different implements of random number generators 
and different random seeds resulting in the generation of completely different sample 
sets. The discrepancy between different MCS experiments is more pronounced with 
small sample size, and could be reduced, but not completely eliminated, by an 
increase, possibly substantial, in the sample size. Therefore, HSS not only leads 
to more accurate estimations than MCS for a given number of samples but also 
eliminates the randomness from the estimation. 

5.2 Mean and Variances 

The reliability bounds and the robustness metrics introduced above depend exclu- 
sively on means and variances. In this paper, Sampling and the First and Second 
Moment Second Order approximations are used for their estimation. These methods 
are briefly introduced next. 
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Figure 3. Points and samples via MCS and HSS. 


5.2.1 Sampling 

Unbiased estimators for the mean and variance of the random variable x are 


1 . . 

E[x] ~ ^ Xi ( 28 ) 

i = 1 

v M w ^i>- E M ) 2 < 29 ) 


where Xi = x(pi) is the ith sample. The p* sample of / P (p) can be generated by 
any sampling technique. 


5.2.2 First and Second Moments of a Second Order Taylor approxima- 
tion (FSMSO) 

These approximations result from calculating the first and second moment of a 
second order Taylor expansion of x(p) about E[p] . Let p E M m , and p, be the ith 
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component of p. For an uncorrelated / P (p), the resulting approximations are 


E[x] 

V[x] 




+£ 


1 

4 



(F[p ( i)] — V[ P (i)] 2 ) 



d 2 x 

0P(i)0P(j) 


) V[p (i) ]V[p u) ] 


T[P(i)] 


(30) 

(31) 


where the functions and derivatives are evaluated at E[p], T[-] is the third central 
moment operator and F[-] is the fourth. 


5.3 Failure probabilities 

In general, reliability metrics cannot be evaluated exactly since they involve the eval- 
uation of complicated integrals, usually multi-dimensional, over complex domains. 
The estimation of failure probabilities, which are basic components of the reliability 
metrics, can be done using sampling or asymptotic approximations. They are briefly 
introduced next. 


5.3.1 Sampling 

The estimation of failure probabilities via sampling is given by 

p f ~ y' 1(. x i g P) 


(32) 


where T(-) is a binary indicator function that gives one if its argument is true 
and zero otherwise. MCS and HSS can be used to generate the required samples 
•El 7 ^2? •••? 


5.3.2 First- Order- Reliability-Method (FORM) 

FORM [12] uses an asymptotic approximation for the estimation of failure probabil- 
ities. In the process, p is transformed into the standard normal uncorrelated space 
q. If p = T(q) = 1 (F q (q)), Equation (6) is equivalent to 

Pf = f ,/q(q)dq 

•O(T(q))<0 

FORM approximates the domain g (T(q)) < 0, by a half-space fitted to the true 
domain at the point of maximum probability density. This approximation leads to 

Pf ~ $(-||q*ll) (33) 
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where q*, called the Most Probable Point (MPP) of failure, is given by the solution 
to the constrained optimization problem q* = argmin{||q|| : g(T(q)) = 0}. In this 
expression, $ refers to the CDF of a standard normal random variable. Notice 
that the rotational symmetry of / q (q) leads to the one-dinrensional approximation 
in Equation (33). When a point satisfying g(T(q)) = 0 does not exist, the MPP 
does not exist and the probability of failure is zero or one. Even though FORM is 
extensively used in structural engineering, its application to controls has been limited 
to stability analysis [11]. Non-smooth limit state functions happen in examples of 
interest. For instance, when the value of A, the maximum real part of the eigenvalues, 
is equal to the real part of more than one of the eigenvalues (counting multiplicity) . 
This consideration shall be taken into account when setting the algorithm to search 
for the MPP. Besides, if the MPP falls at a derivative discontinuity on the limit state 
surface, the approximation given by Equation (33) is not as good as first order. 

5.3.3 Hybrid Approach 

Sampling based techniques can readily be used to estimate probabilities of failure 
using Equation (32). However, high computational demands in the evaluation of 
Xi = x(pi) can preclude their practicality especially when Pf ~ 0 (or Pf ~ 1). 
Examples of this can be easily found 3 . On the other hand, methods based on 
asymptotic approximations, such as FORM, provide good approximations when Pf 
is small. This is clear since for large values of q, / q ( q) decreases exponentially in 
||q|| 2 , so if q* is large enough, most of the failure probability comes from the part 
of the failure event in the immediate neighborhood of q*. On the other hand, for 
failure probabilities away from zero and one, the slow decrease in F q (q) near the 
MPP and the geometrical difference between the true limit state function and its 
linear approximation contribute a bigger error to the FORM approximation. 

In this paper, a hybrid approach which combines HSS and FORM is used to 
estimate probabilities of failure. In order to identify the numerical tool that best 
suits the task at hand, a coarse and computationally-efficient estimation of Pf is first 
generated using HSS. Such estimation is then compared with a reference, namely the 
reference failure probability P re f, a user-defined value set in advance. The compari- 
son between the coarse estimate and P re f is used to determine which of FORM or 
HSS will be used to generate a new estimation, presumably more accurate. Assume 
that two sets of Hammersley samples of / p ( p) are available. One set has n\ samples 
and the other one has ri 2 samples, where 712 3> n±. For a given failure domain T 
and a user defined reference failure probability P re f, proceed as follows 

1. A 1 Assessment : estimate the Pf using Equation (32) and the set of n± samples. 

2. A2 Assessment: recalculate Pf as follows. If the estimated value is greater 
than P re f use Equation (32) with the set of n2 samples. If the estimated value 
is less than P re f use FORM. 

* s Wang and Stengel [6, 7] make the approach computationally viable by using a single random 
variable to model 28 uncertain parameters. The same authors [14] require 25000 samples to deter- 
mine a sufficiently small 95% confidence interval. 
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The refinement performed in A2 might not always be necessary. Situations in which 
this is the case are as follows. Since reliability metrics for random processes are 
heavily dependent on the larger values of the probabilities of failure that compose 
them, (see the bottom plot of Figure 2 ), refining the estimation of the small failure 
probabilities is inconsequential. Furthermore, if the reliability metrics are used to 
calculate the cost function of an optimization problem, more accurate estimations 
are not needed when the assessment resulting from using the coarse estimate of Step 
1 indicates a poor design, e.g. Fa(0) 3> 0. 


6 Control Synthesis 

6.1 Reliability-based 

The formulation of the control design problem from a reliability perspective is as 
follows. For a given plant model, compensator structure, uncertainty model and a 
set of design requirements prescribed via inequality constraints; one would like to 
find the compensator parameters for which the resulting probability of violating the 
design requirements is minimized. Notice that this refers to the excursion of the 
outcomes into the failure domains. 

6.1.1 Robustness Considerations 

The reliability metrics in Equations ( 7 - 10 ) are usually applied using a fixed failure 
set T. In this form, a reliability analysis cannot assess the system’s performance in 
the regions where the design requirements are satisfied, i.e. , the intersection of the 
admissible domains associated with all the design requirements. Since the portion 
of the random outcome lying in the admissible domain A might end up being sub- 
stantially larger than the portion lying in the failure domain T .. a reliability-based 
approach with fixed failure boundaries does not have control over the bulk portion 
of the PDF, which is the portion that dictates the most likely performance. 

We now introduce the concept of a shapable failure set T with an example. Let 
,x(k) be the stationary RMS value of an error signal. One would like to find k such 
that x is as close as possible to zero. Uncertainty in the plant makes x a random 
variable. Let x be the failure boundary associated with the design requirement 
x < x, i.e., a fixed failure set is T = {x | x E [x, 00)}. The minimization of 
r x (x) leads to reliability optimal compensators. Suppose there exist multiple designs 
leading to r x (x ) = 0 . These designs however differ in how well the resulting PDF 
of x spreads over the admissible domain A = {x \ x E [0, Zc) } . The concentration 
of f x (x) about zero is an indicator of the robust performance. Say, ki leads to 
f x (x/ 2 ) = 0 which implies that r x (x ) = 0 while U.2 leads to f x (x) = 0 but f x (x/ 2) > 
0 . Since neither of these two designs violate the design requirement x > x, a 
reliability analysis cannot establish that the compensator with parameters ki has a 
better robust performance than the one which uses k2. By minimizing the reliability 
metrics and simultaneously shrinking the admissible domain, the uncertain system 
performance can be concentrated in a more desirable region. This is attained by 
parameterizing the failure boundaries of J- as well as a penalizing function 7^ with 
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Figure 4. Robust performance concepts and shapable failure domains 


an additional design variable, namely e. The basic idea is to solve an optimization 
problem for the design variable d = [k, e] such that a twofold objective is pursued: 
a reliability metric is minimized while the size of A is reduced. This implies that 
the failure and the admissible domains are now dependent on e, i.e. , -4(e) and 
.T(e). For the RMS example above, the minimization of J = r x (e) + ^y x where 
= e, d = [k, e] and e E [0,x], leads to the desired solution. This setting implies 
lF(e) = {x | x E [e, oo)} and -4(e) = {x \ x E [0,e)}. See Figure 4. Notice that the 
value of J for ki is less than the one for k 2 if e E \x/2,x). 

In general, we will refer to the augmented reliability metric as the sum of a reli- 
ability metric from Section 3 and the penalizing term introduce above. Augmented 
reliability metrics for the random variable x and the random process x(h) take the 
form r x (x(e), x(e)) + 7 x (e) and r x ^}A{x(h,e),x(h,e)) + 'Y x (h)( e ) The penalizing func- 
tions 7 a; (e) and 7 x (h){ e ) must be proportional to the size of the admissible domain 
-4(e). In addition, they must be built such that the minimization of the augmented 
metric does not lead to unacceptable solutions, e.g., r x = 1 and 7 x = 0. If r x < e is 
required, use a monotonically increasing function satisfying 7 E [0 , e] . In the RMS 
example above, for which x < x is the design requirement, this is attained by min- 
imizing an augmented reliability metric with 7 r(e) = ee/x for e E [0, x] , x(e) = 0 
and x(e) = e. 

6.2 Robustness-based 

The formulation of the control design problem from a robustness perspective is as 
follows. For a given plant model, compensator structure and uncertainty model; 
one would like to find the compensator parameters for which the resulting random 
outcome is as concentrated as possible around a target deterministic behavior. While 
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the deterministic target for a random variable is a scalar, the target for a random 
process is a function, e.g., u(t) = 0 is the target function for the control. 

6.3 Multi-objective Optimization 

In order to satisfy multiple design requirements, all the corresponding reliability 
metrics must be driven to zero. Since some of these objectives are conflicting, a 
compromise among them might be required. The technique to be used reduces the 
multi-objective problem into a single objective problem given by a weighted sum 
of the reliability metrics. The solution to this problem will converge to a Pareto 
optimal point, which is a design for which no improvement on a single objective can 
be attained without making another worse. For each set of weights, the optimiza- 
tion may converge to a different Pareto point. Other techniques, such as the goal 
attainment method could be used instead. 

6.4 Synthesis Procedure 

A step-by-step procedure for control synthesis is presented next. 

1. Determine the plant model and the control structure. First principles and 
classical deterministic approaches to compensator design can be used. Identify 
the set of parameters that have a strong impact on the plant model. Use 
sensitivity information and engineering judgment to select the set of uncertain 
parameters p. At this stage, the parametric plant model, i.e., G*(p), and the 
control structure, i.e., K{ k), must be fully determined. 

2. Generate the probabilistic parameter model / P (p). Use engineering judgment 
and experimental data if available. 

3. Use Equation (26) to generate the sample sets of / P (p) for both n\ and ri 2 - 

4. Cast the design requirements, either reliability or robustness based, in terms 
of the metrics introduced above. Use Equations (7,10) for the reliability met- 
rics and Equations (24,25) for the robustness metrics. Recall that while each 
reliability requirement requires setting a failure domain T ’, each robustness re- 
quirement requires setting a target behavior. Use these metrics to compose the 
cost vector c, which is the vector of objectives for multi-objective optimization. 
Recall that robustness considerations can be alternatively considered within a 
reliability formulation, as explained in Section 6.1.1. 

5. Let d be the design variable. For robustness-based metrics and reliability- 
based metrics with fixed failure domains, d = k. Reliability-metrics with 
shapable failure domains lead to d = [k T ,e T ] T . Build a penalizing function 
7 (e) for these terms following the guidelines provided in Section 6.1.1. Update 
the components of the cost vector c by adding the penalizing functions and 
parameterizing the failure boundaries. 
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6. Solve the single objective optimization problem 

d* = argrnin { c T w } 


( 34 ) 


where w is composed of non-negative weights. Each cost function evaluation 
used in the search for the optimal design d* requires a probabilistic analysis. 
This analysis is done by calculating the metrics contained in c. While the 
hybrid approach is used for estimating the reliability metrics, either HSS or 
FSMSO are used for the robustness metrics. This task requires forming the 
closed-loop Equations (4-5) and performing typical control studies such as 
finding closed loop poles, time responses and Bode plots. 

During optimization, the following procedure is suggested in order to focus 
most of the computational effort toward the assessment of better designs. For 
a reliability metric, first perform the A\ assessment of Section 5.3.3 for n\ 
samples of p and e\ h-samples. This implies that only the first step of the 
hybrid approach is applied to all reliability metrics. If A\ shows that d is a 
good design, perform the refined assessment A 2 . A 2 is carried out by using 
n 2 samples of p, e 2 h-samples and a user defined value for P re f- The rational 
of this was presented in Section 5.3.3. For a robustness metric or a reliability 
bound, HSS or FSMSO could be used. In the former case, a first assessment 
can be done with the parameters of A\ and a second one with the parameters 
of A 2 . This two-fold analysis procedure is used in the examples. 

The cost vector c can be formed by combinations of reliability-based metrics, relia- 
bility bounds and/or robustness-based metrics. The corresponding implications are 
explored in the examples. 

6.5 Optimization under Uncertainty 

Let J(p,d) = c r w denote the cost function of the optimization problem. Due to 
the nature of the reliability metrics in c, the cost function might not only have 
plateaus, i.e., there could exist a design d and a non-zero perturbation 6 such that 
J(p,d) = J(p,d + 5), but might also have a discontinuous gradient. 

The use of sampling in the estimation of probabilities makes the cost function 
piecewise constant. Let J e (p,d) be an estimation of the actual cost J(p,d). For 
any design d and regardless of the number of samples, there always exists a pertur- 
bation 5 such that J e (p, d) = J e (p, d + 5). This situation is aggravated, i.e., bigger 
perturbations can be found, when a smaller number of samples is used or when Pf 
is close to zero or one. 

Mean and variances can be estimated via sampling or via FSMSO. Their es- 
timation via sampling does not exhibit the numerical problems mentioned above. 
Among the sampling techniques, it was found that HSS suites very well the need 
for an accurate and efficient estimation [18]. Estimation via FSMSO is considerably 
faster than sampling but less accurate in general. This is so, since the Taylor ex- 
pansion might become a poor approximation of the actual function when / P (p) is 
not concentrated about E[p], 
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The nature of J e must be taken into account when selecting a numerical method 
for optimization. In the examples to follow, the resulting non-convex piecewise con- 
stant optimization problem is first solved using Genetic Algorithms (GA) for a fixed 
number of generations. Since GA is based on a random search, the hybrid approach 
is particularly convenient. After the fixed number of generations is reached, the 
GA solution is refined using the Nelder Mead Simplex algorithm, which is a local 
non-gradient based search method. 


7 Numerical Examples 

The synthesis procedure of Section 6.4 is applied herein to two examples. A textbook 
satellite attitude control problem is considered first. Then, the active control of a 
flexible beam subject to disturbances is considered. The following parameter values 
are assumed in both problems. If p € M m , the coarse assessment Ai is performed 
using m = 75 m p-samples and ei = 90 h— samples. For the finer assessment A 2 , 
n 2 = 500m p-samples, e 2 = 180 h — samples and P re f = 0.01 are assumed. For the 
sake of comparison, the examples also present the analysis of deterministic versions 
of the problems for which E[p] is used as parameter. Such problems and their 
solutions will be referred to as the nominal ones. 


7.1 Satellite Attitude Control 


Accurate satellite pointing in the presence of large thermal gradients and mass losses 
for uncertain initial conditions is desired. A simple rotational model of two bodies 
connected with a flexible boom leads to 


J\ 6 \ + b( 6 i — 62 ) + k(9i — 62 ) = u 
J 2 O 2 + b(0 2 — 9\) + k{9 2 — 61 ) = 0 


where 6 \ and 62 are the deflection angles, J\ and J 2 are moments of inertia, k is 
the equivalent stiffness, b = ayjk / 10 is the equivalent damping coefficient and u is 
the applied torque. The variable a is used to model the changes in damping caused 
by thermal variations. Assume that J 2 = 0.1 since mass losses only affect J\. The 
non-collocated sensor-actuator pair resulting from using y = 62 leads to the SISO 
system 


G( P ) 


k + bs 

J 1 J 2 S 4 + b(J\ + J 2 )s 3 + (Ji + T 2 ) fes 2 


(35) 


Variations in the operating conditions and ignorance on the initial conditions are 
modeled using p = [J\,a, k, 9i t o, 0i,o, 02, 0 , # 2 ,o] T , where 6b , 0 = 0i(O) and 0 2 , 0 = # 2 ( 0 ). 
The following output-feedback control structure is assumed 

k\ + k2S + kgS 2 + /’4-S 3 

kg + kes + kjs 2 + kgs 3 




The joint PDF that describes the uncertainty in p is given by the independent 
random variables listed in Table 1, where U and B refer to Uniform and Beta 
distributions. 
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Table 1. Uncertainty Model. 


Jl 

Aji = [0.8, 1] 

fM) = 17(0.8,1) 

a 

A a = [0.03,0.2] 

f a (a) = B( 0.3, 0.2) 

k 

A k = [0.09,0.4] 

f k (k) = B( 5,5) 

0i,o 

A 01iO = Hr/2, tt/2] 

fg l o (0 1 ,o) = B( 5.2, 5.2) 

0i,o 

^ei.o = [-15,15] 

4, 0 (M = 5(2.5, 2.5) 

@2,0 

A <b,0 = [-ttAtt/ 2 ] 

fe 2 , 0 (@ 2 ,o) = 5(5.2,5.2) 

@2,0 

A feo = [— lS.lS] 

4, 0 (M = 5(2.5, 2.5) 


7.1.1 Nominal Compensator 

A baseline compensator for the nominal plant is designed by standard pole place- 
ment techniques such that large stability margins are attained. For fairness sake, the 
baseline compensator’s gains where selected such that the resulting closed-loop sys- 
tem is robustly stable for the uncertainty model of Table 1. The resulting gains are 
ki = 10 6 [0.0108, -0.3271,0.1192, 0.0092,1.8835,2.1305, 2.2276, 0.9308] T . The anal- 
ysis of the nominal compensator using / p ( p) indicates that the closed-loop system 
is robustly stable as intended, i.e. r\(0) = 0, but the time responses are unsatisfac- 
tory. The CDF of A as well as the time evolutions of the output and the control 
signals are shown in Figures 5-7. The sudden variation in the slope of the CDF 
of A in Figure 5 is the result of a change in the closed- loop pole that determines 
A. The considerable disparity between A(E[p]) and E[A(p)] shows that the nominal 
problem is not a meaningful representative of the probabilistic behavior. Figures 6 
and 7 show the time evolution of the random signals by indicating the 1, 10, 20, 
30, 40, 50, 60, 70, 80, 90 and 99 percentiles. In Figures 6 and 7, the percentiles 
and the nominal functions are shown. Dotted lines are used to indicate the failure 
boundaries introduced in the next section. It is interesting to see how the PDFs 
expand, e.g. Figure 7 at 2.5 and 8 seconds, and contract, e.g. Figure 7 at 4 and 
16 seconds, in a oscillatory manner. This information can be used to determine the 
time periods when the effects of uncertainty are more noticeable. 

7.1.2 Reliability-based compensator 

Design requirements for a reliability formulation are as follows. Performance require- 
ments on the closed-loop stability, settling time, over-shoot, control usage and on the 
magnitude of the loop transfer function lead to c = [r^(A) , r y ^ ( y{t ) , y(t)) , r u ^ ( u , u) , 
Uj( a ,)(g(w))] T , where q(ui) = \GK\ is the loop gain. The failure bound- 
aries to be used are A = 0; y(t) = —1.257i(t) + 2.27i(t — 70) for t G [0, 80], where TL is 
the Heaviside function; y(t) = 1.2577(f) — 0.277(f — 70) for t G [0,80]; u(t) = —0.5 for 
t G [0,25]; u(t) = 0.5 for t G [0,25]; q(co) = 0.75 / c o for u G [10 — 6 , 0.2] and q{u) = 1 
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Figure 5. CDF of A for the nominal compensator. 


- Percentiles 


y(t) at E[p] 




Figure 6. y(t) for the nominal compensator. A zoom is shown below. 
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Figure 7. u(t) for the nominal compensator. 


for io G [1, 10 2 ]. 

The synthesis algorithm for a fixed set of failure boundaries leads to kg = 
10 6 [0.0405, 0.1267,0.2422,0.0320, 0.5244, 1.0057, 1.2263, 0.6560] T and c = [3.13 x 
10 _4 ,6.51 x 10 -4 , 0.0037,6.65 x 10 _4 ,0] T for which the weighting vector w = 
[500, 1, 1, 1, l] r was used. This weighting vector was selected in order to empha- 
size stability. A probabilistic analysis of this compensator leads to Figures 8-11, 
where the CDF of the dominant closed-loop poles, the output, the control and 
the loop gain processes are displayed. Better robust stability characteristics, i.e., 
smaller values of r\(0), in this compensator are attained by increasing w^p Recall 
that reaching zero probability of instability might be unfeasible. From Figure 10 
we see that for all possible parameter values and initial conditions the process u(t ) 
stays within the ±0.5 range with more than 0.8 probability after 6 seconds. Fig- 
ure 11 shows that uncertainty mostly affects the damping and the location of the 
resonant frequency. While all the percentiles for u> < 0.1 are indistinguishable from 
the nominal curve, percentiles in the ui G [0.1, 1] are not shown. In contrast to the 
mid frequency range, where the effects of uncertainty are large, there is no violation 
of the low frequency design requirement. The CDFs of A for the nominal and the 
reliability based compensator are superimposed in Figure 12. Despite the increased 
variability of the dominant closed-loop poles of the reliability-based compensator 
(the support of A is about two and a half times larger than the one for the nomi- 
nal compensator), the system is robustly unstable with 3.13 x 10~ 4 probability. A 
substantial improvement in the performance is achieved by trading off a very small 
margin of the probability of stability. This improvement can be seen after comparing 
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Figures 6-7 with 9-10. Overall, the performance resulting from k 2 is substantially 
better than the one resulting from ki. 

During optimization, 157 random variables were used to evaluate c for the coarse 
assessment A\. This task took 23.6 seconds when performed on a Pentium III 1795 
MHz with 512 MB of RAM. Notice that the CPU time associated with A 2 depends 
on the initial conditions used to find the MPPs. For this assessment, HSS was used 
for 628 random variables and the hybrid approach was used for robust stability. 
This task took 102 seconds. 


E[X(p)] *<E[p]) 



Figure 8. CDF of A for the reliability-based optimal compensator. 


7.1.3 Using Mean and Variance based metrics 

If only means and variances are used, reliability and robustness metrics can be com- 
bined to form c = [b A (0), r y ^( 1), r u(t )(0), 6^(1), 6 g ( t y ) (0.75/u;)] T . Ranges provided 
before will be used with the exception of u; G [5, 10 2 ] for b. All the components of 
the cost vector require of the estimation of means and variances, task that can be 
performed via sampling or the FSMSO method. 

HSS leads to k 3 = 10 5 [1.1207, 0.2307, 1.5235, 0.1713, 0.001, 0.0692, 1.2294, 1.069] T 
and c = [8.6 x 10~ 5 , 0.0057, 0.0129, 1.94 x 10~ 4 , 6.189x 10~ 5 ] r for which the weighting 
vector w = [500, 1,4, 1, 1] T was used. A probabilistic analysis of this compensator 
leads to the results shown in Figures 13, 14 and 15. Figure 13 shows that the opti- 
mal compensator makes A insensitive to uncertainty, i.e. the CDF resembles a step 
function which would result if the system is deterministic. Notice that the bound on 
the stability condition tends to reduce the variability of A. This artificially reduces 
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Figure 9. y{t) for the reliability-based optimal compensator. A zoom is shown 
below. 


the design space, i.e. there might exist solutions with a better performance in which 
large variances in A occur. This is a consequence of the form of the bound. The 
comparison between the time responses for the reliability compensator and Figures 
14 and 15 show substantial differences between the two solutions. While the control 
for the reliability based compensator intends to keep the process within the strip 
|n| < 0.5, the mean and variance based solution intends to concentrate the random 
process about the target function u = 0. It can be seen that this is achieved by 
inducing oscillations about the target function. The reader should also notice that 
excursions beyond the failure boundaries in a reliability formulation, e.g. y(t) in 
Figure 9 for the first 10 seconds, are not penalized according to the severity of the 
violation. This is in sharp contrast with the mean and variance formulation. Since 
we are using the same number of samples as before, the savings in CPU time when 
comparing this formulation with the one of the previous section result from not us- 
ing FORM. Hence, the CPU time for A\ is about the same while the CPU time for 
A 2 is 77s. For a given number of samples, the estimation of means and variances is 
in general more accurate than the estimation of failure probabilities. Small sample 
sets however, result in larger estimation errors of the moments, which could lead to 
the selection of the wrong conditions in Equations (19-20) and (22-23). 

Next, the FSMSO method is used for estimating the cost vector. The first and 
second order derivatives for all the metrics of interest, i.e., closed-loop poles, output, 
control and loop transfer function, were derived analytically. Some of the required 
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Figure 10. u(t) for the reliability-based optimal compensator. 



Figure 11. Bode plot of the loop gain for the reliability-based optimal compensator. 
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Figure 12. CDFs of A. 
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Figure 13. CDF of A for the mean and variance based compensator. 
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Figure 14. y(i) for the mean and variance based compensator. A zoom is shown 
below. 


expressions are as follows 
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where A, B, C, and D were given in Equations (4) and (5), x = [x^,x T ] T , <9[-] 
indicates a derivative with respect to p(j) excluding initial conditions, and 
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Figure 15. u(t) for the mean and variance base compensator 


B=[B C BD C ] T , dB = [ 0 dBB c ] T , d 2 B = [ 0 <9 2 BD C ] T , 

C = [ 0 C ] , dC = [ 0 dC ] , d 2 C = [ 0 d 2 C } 

In these expressions, the subscript c refers to the state space representation of the 
compensator I\ while the matrices with no subscript refer to the open-loop plant. 
The time evolution of the sensitivities is calculated by solving the state space model 
in Equations (37-38). Sensitivities with respect to initial conditions can be ana- 
lytically calculated via the matrix exponential. Such developments as well as the 
ones for the other sensitivities are omitted due to space limitations. In general, the 
analyses based on FSMSO were inaccurate when compared with sampling. Even 
though the CPU time per analysis was reduced to 2 seconds, large errors in the 
estimation precluded its use for synthesis. Table 2 shows a comparison between 
the assessments resulting from HSS and FSMSO as the variance of the uncertain 
parameters is increased. The same PDFs for the input uncertainty of Table 1 are 
used while the means are kept constant. The right most column shows the aver- 
age relative error in the components of c. For the Beta distributions, the increase 
of the variance was attained by enlarging the support. The terms resulting from 
the crossed derivatives were neglected. It can be observed that the accuracy of the 
FSMSO estimation rapidly decreases as / P (p) is less concentrated about its mean. 
This trait is obvious since large excursions from E[p] might considerably degrade 
the accuracy of the Taylor approximation. 
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Table 2. Comparison of HSS and FSMSO for k 3 . 


V[P(,:)] for 
i = 1, ... 7 

Cost vector 

c 

Average 

error 

1 x 10 -4 

c H ss = 10~ 3 [0.000823, 0.0096, 0.5880, 0.0058, 0.0013] 
cfsmso = 10” 3 [0.000827, 0.0095, 0.5880, 0.0058, 0.0001] 

0 . 3 % 

1 x 10~ 3 

chss = 10“ 3 [0.00824, 0.0110, 0.6560, 0.0580, 0.0136] 
cfsmso = 10~ 3 [0.00823, 0.0108, 0.7720, 0.0580, 0.0133] 

4 % 

3 x 10~ 3 

c H ss = 10~ 3 [5.1000, 0.0505, 0.8800, 0.1745, 0.0446] 
CFSMSO = 10" 3 [0.0247, 0.0194, 1.9160, 0.1740, 0.0408] 

57 % 

4 x 10~ 3 

chss = 10“ 3 [32.400, 0.0013, 1.0880,0.2330,0.0626] 
cfsmso = 10" 3 [0.0329, 0.0267, 2.8680, 0.2320, 0.0548] 

454 % 


To show the reduction in the design space caused by using b\ (0), we have solved 
this problem using f A (0) instead. This practice results in the compensator k.j = 
10 5 [0.9082, 0.3590, 1.1773, 0.2372, 0.001, 0.0616, 1.0432, 1.196] T for which c = [1.86 x 
10~ 5 , 0.0052, 0.0121, 1.845 x 10~ 4 ,6.10 x 10“ 5 ] T . Figure 16 shows that A a is about 
10 times larger than the one in Figure 13. Actually, V[A] is about 2276 times larger. 
This fact heavily penalizes the compensator via b\( 0) in spite of its actual improved 
performance. The reader should also notice that there is an offset of about 10% 
between E[A(p)] and A(E[p]). Performance improvements in all the metrics were 
attained. Figures 17 and 18 show the corresponding output and control, from which 
considerable improvements are seen when compared with Figures 14 and 15. 

7.2 Disturbance Rejection for a Flexible Beam 

The second example will focus on the disturbance rejection for a flexible beam test 
article with both physical and modal parameter uncertainties. The system, displayed 
in Figure 19, consists of a flexible thin aluminum blade, one-meter long, attached at 
its base to a hub motor. The hub motor is the control actuator for the system. At 
the tip of the beam, there is a reaction wheel that serves as a disturbance generator. 
The test article has nine sensors that may be used in any combination for either 
feedback or performance output monitoring. The finite element method is used to 
model this system by utilizing Euler-Bernoulli planar beam elements. A complete 
description of the flexible beam test article [20] is available. 

For this paper a SISO system is studied. The input u is the hub motor torque 
and the measured output y is the tip velocity. The tip reaction wheel disturbance 
is modeled by passing a Gaussian white noise process through a second-order linear 
low-pass filter, with parameters £/ = 0.8 and toj = 2007T rad/s. The first five modes 
of the elastic structure are used to build a state space realization of the plant. This, 
in addition to the disturbance model leads to an open-loop system where x £ M 12 , 
u £ M and y £ M. The uncertain parameters are the Young’s Modulus E (Pa), 
the density p (Kg/rn 3 ) and the damping ratios of the retained vibration modes £j, 
i = 1,2, 3, 4, 5. This set leads to p = [E, p, £i, £ 2 , £ 3 , £ 4 , £s] T , whose components are 
assumed independent. The corresponding PDFs are given in Table 3. The mean 
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Figure 16. CDF of A for a compensator with parameters k 4 . 
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Figure 17. y(t) for a compensator with parameters k 4 . A zoom is shown below. 
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Figure 18. u(t) for a compensator with parameters k 4 



Figure 19. Flexible beam test article. 
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Table 3. Uncertainty Model. 


E 

A e = 10 lo [5.226, 7.839] 

f E (E) = B( 5,5) 

P 

A p = [2280, 3420] 

fM = B{ 3,3) 

6 

A ?1 = [0.08,0.12] 

fti(Zi) = B( 2,2) 

6 

A 6 = [0.0252,0.0378] 

/&(&) = £( 2,2) 

£3 

A fo = [0.02,0.03] 

/&(&) = £(2,2) 

U 

A ?4 = [0.0304,0.0456] 

/&(&) = £(2,2) 

£5 

A ?5 = [0.02,0.03] 

/&(&) = £(2,2) 


values of the uncertain parameters are set to coincide with parameters in the finite 
element model leading to good matches with experimental data. The supports of 
the distributions were set according to reasonable ranges of variation. The shapes 
of the PDFs were chosen arbitrarily. Performance requirements on stability and the 
output RMS are considered. Full-state feedback with a full-order observer determine 
the control structure. 

7.2.1 Nominal Compensator 

As before, a baseline compensator for the nominal plant is designed such that the 
RMS value is minimized. The resulting compensator, with parameters di, leads 
to Urms = 0.011 m/s. The propagation of the uncertainty prescribed in Table 3 
through the closed loop-system show that this compensator is robustly unstable 
with Tx(0) = 0.235. 

7.2.2 Reliability-based Compensator 

For this example, a shapable failure domain for the RMS requirement is assumed. 
This leads to the cost vector c = [r : A(0),r 2/rms (e) + 7j /rms ] T , where e 6 [0,0.05] and 
7 yrm S = e. The selected control structure makes the feedback gain G, the observer 
gain L and the RMS failure boundary e, the design variables. Recall that the 
separation principle does not hold. The resulting closed-loop dynamics is given by 
Equations (4) and (5). Notice that although the observer is deterministic, all the 
closed-loop poles are random. 

The synthesis approach with w = [20, 1] T leads to a compensator with pa- 
rameters d 2 , for which f\(0) = 0, e = 0.0139 m/s, r yrms (e) = 3.6 x 10 -3 and 
c = [0, 3.6 x 10' 3 ] T . The probabilistic analysis of this compensator leads to the re- 
sults shown in Figures 20-21. Figure 20 shows that the random variable y rms for d 2 
is moved toward zero from y rms = 0.05, by virtue of the non-fixed failure boundary. 
Figure 21 shows Bode magnitude plots of the disturbance to output transfer func- 
tion, namely T zy . Notice that differences in the mid-frequency range, i.e. u ~ 1, of 
the diagram have a bigger impact on the RMS value. In addition, considerable vari- 
ability in the closed-loop Bode magnitude plot as well as a significant reduction in 
the damping of the first mode are attained near 20 rad/s. It is interesting to notice 
that even though d 2 leads to a robustly stable closed-loop system in Equation (4), 
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the full-state feedback subsystem Agi and the full-order observer subsystem A- 2,2 
have a non-zero probability of instability. This indicates that the application of the 
Separation Principle before accounting for uncertainty in the model artificially re- 
duces the design space. In other words, a design that accommodates for uncertainty 
based on the full system dynamics, i.e. Equation (4), may lead to robustly stable 
solutions in a design region which would have been rejected for reasons of instability 
if the Separation Principle had been applied before searching for a robust regulator 
with full-state feedback and a robust observer. 



Figure 20. PDFs for the reliability based and the robustness based compensators. 


7.2.3 Mixed Compensator 

Lets take c = [r^(0), T yrrn3 (0)] T , where the hybrid approach will be used for the 
stability metric and HSS for the RMS metric. The synthesis algorithm leads to 
d 3 , for which c = [0,1.34 x 10~ 4 ] r . This compensator leads to the dashed line in 
Figure 20. Comparing both solutions, we see that while the PDF corresponding 
to the reliability-based compensator has less probability of exceeding the boundary 
value of e = 0.0139 nr/s, the robustness-based compensator leads to a PDF which is 
much more concentrated toward the ideal value of zero. This clearly shows that the 
conceptual differences between the two formulations. Since there is no conservatism 
in the selection of the nominal plant, i.e., G(E[p]) is not the most difficult plant to 
control, yrms = 0.011 m/s does not necessarily bound the supports of the PDFs. 
This can be observed in Figure 20, where the support corresponding to d 3 contains 
yrms = 0.011 m/s. 
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Figure 21. Bode diagrams of T zy for d 2 - 


7.2.4 Mean and Variance based Compensator 


For this case we assume c = [6^(0), %, ms (0)] r , where the FSMSO method will be 
used to estimate both indexes. In contrast to the previous example, the accuracy of 
the results made the method suitable for synthesis. Analytical sensitivities were used 
in the moments approximations. For instance, the sensitivities of the closed-loop 
poles that determine A are given by 


Av i = 

dsj = zJdAvj 


d 2 Sn = z i 


d 2 A + (dA - dsjlXsjl - A) <9A + dA( Sj I - A) (<9A - d Sj I 


where z j is the jth eigenvector of A 1 , the right and left eigenvalues are normalized, 
i.e. z Jvj = 1, and and non-repeated poles are assumed. Derivatives of the output 
covariance are given by 


dy rms = jdiag C<9QC T + 25CQC 5 


1 1/2 


d 1 


y rms = |diag C9 2 QC T + 4:dCdQC T + 2d 2 CQC T + 2dCdQdC T } 


1/2 


where the derivatives of the state covariance are given by the solution to the set of 
Lyapunov equations 


AdQ + dQA T + dAQ + QdA T + dBSB T + BSdB T = 0 
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A3 2 Q + d 2 QA T + 2dAdQ + 2dQdA T + d 2 AdQ+ 
dQd 2 A T + 2dBSdB T + d 2 BSB T + BSd 2 B T = 0 

The synthesis algorithm leads to a compensator with parameters d 4 , for which 
c = [5.33 x 10~ 6 , 1.33 x 10 -4 ] T . The resulting PDF for the RMS is indistinguishable 
from the one shown in Figure 20 even though the compensators are different. The 
CDFs of A for d 3 and d 4 are superimposed in Figure 22. As before, the tendency of 
the formulation of making A as deterministic as possible is apparent. Notice however, 
that the CPU time required for the analysis of d 3 is 421 seconds while the one for d 4 
is 1.22 seconds. On the other hand, the analysis of the compensator with parameters 
d 4 via HSS takes 64 seconds. This exemplifies substantial savings in CPU time which 
result from using the FSMSO method. In this example, those savings justify the 
labor required to compute analytical derivatives. The reader should recall, however, 
that the same method led to inaccurate results for the satellite example. 
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Figure 22. CDFs of A for d 3 and d 4 . 
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8 Conclusions 

This paper proposes a control synthesis methodology for systems with probabilis- 
tic uncertainty. Synthesis is performed by solving a multi-objective optimization 
problem which combines requirements of stability and performance in time- and 
frequency-domains. In this study, reliability- and robustness- based formulations 
are proposed and several numerical methods for estimation are examined. In a re- 
liability formulation, the probability of violating design requirements is minimized 
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while admissible domains are contracted toward regions with an improved perfor- 
mance. In a robustness-based formulation, a metric that measures the concentration 
of the random variable/process about a target scalar /function is minimized. These 
two formulations lead to compensators with distinctive characteristics. In addition, 
metrics that bound the reliability metrics proposed, whose estimation only requires 
of means and variances, are also derived and used for control design. 

Some of the fundamental differences between the proposed strategy and conven- 
tional robust control methods are: (i) unnecessary conservatism is eliminated since 
there is not need for convex supports/sets, (ii) the most likely plants are favored 
during synthesis allowing for probabilistic robust optimality, (iii) the tradeoff be- 
tween robust stability and robust performance can be explored numerically, (iv) the 
uncertainty set, which could be unbounded, is closely related to parameters with 
clear physical meaning, and (v) compensators with improved robust characteristics 
for a given control structure can be designed, e.g., one can search for a PID con- 
troller with best robust characteristics. Examples related to the attitude control of 
a satellite and to control design for disturbance rejection in a flexible beam are used 
to demonstrate and validate the methodology. 
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